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A HYBRID SEGMENTATION AND D-BAR METHOD FOR ELECTRICAL 

IMPEDANCE TOMOGRAPHY 

S. J. HAMILTON, J. M. REYES, S. SILTANEN, AND X. ZHANG 


Abstract. The Regularized D-bar method for Electrical Impedance Tomography provides 
a rigorous mathematical approach for solving the full nonlinear inverse problem directly, 
i.e. without iterations. It is based on a low-pass filtering in the (nonlinear) frequency 
domain. However, the resulting D-bar reconstructions are inherently smoothed leading to a 
loss of edge distinction. In this paper, a novel approach that combines the rigor of the D- 
bar approach with the edge-preserving nature of Total Variation regularization is presented. 
The method also includes a data-driven contrast adjustment technique guided by the key 
functions (CGO solutions) of the D-bar method. The new TV-Enhanced D-har Method 
produces reconstructions with sharper edges and improved contrast while still solving the full 
nonlinear problem. This is achieved by using the TV-induced edges to increase the truncation 
radius of the scattering data in the nonlinear frequency domain thereby increasing the radius 
of the low pass filter. The algorithm is tested on numerically simulated noisy EIT data and 
demonstrates significant improvements in edge preservation and contrast which can be highly 
valuable for absolute EIT imaging. 


1. Introduction 

In Electrical Impedance Tomography (EIT) a conductive body is probed with harmless elec¬ 
trical currents fed into the body through electrodes at the surface, and the resulting voltages 
are measured at the electrodes. The goal is to recover the electrical conductivity distribu¬ 
tion inside the body from these surface electrical measurements. EIT is useful in medical 
imaging, as different tissues have different conductivities, and it allows harmless and painless 
monitoring of patients even over long periods of time. Another application area of EIT is 
non-destructive testing. See [miss] for reviews of EIT and its uses. 

The image reconstruction task of EIT is a nonlinear and severely ill-posed inverse problem. 
Therefore, EIT algorithms need to be regularized to overcome the extreme sensitivity to 
modeling errors and measurement noise. Among EIT algorithms, the so-called D-bar method 
stands out due to its unique capability of dividing the measurement information neatly into 
stable and unstable parts in a (nonlinear) frequency domain. See Figure]^ 

Generally speaking, regularization involves complementing insufficient measurement infor¬ 
mation by a priori knowledge about the conductivity. The D-bar method does this very 
explicitly assuming that the conductivity is twice continuously differentiable, which allows 
replacing the values of the nonlinear Fourier transform by zero in the unstable part of the 
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Figure 1. Left: simulated “heart-and-lungs” phantom conductivity. Middle: 

D-bar reconstruction based on nonlinear low-pass filtering with 0.75% relative 
noise added to FIT voltage data (see Figure [^. Right: reconstruction with 
the proposed hybrid method from the same noisy FIT data. 

frequency domain. Therefore, this low-pass filtering has the side effect that the resulting 
D-bar reconstructions are always smooth, as seen in the middle image of Figure [T] 

In many applications of FIT, including medical imaging, it is important to see boundaries 
between regions of different conductivities. The standard D-bar method practice of inserting 
zeroes for high frequencies is not ideal as crisp boundaries between different conductivity 
regions necessarily contain high frequencies. 

We introduce a novel edge-enhancing method for FIT, built upon the assumption that 
we know a priori that the conductivity is piecewise constant. It builds upon the stable D- 
bar reconstruction and increases the radius of reliable scattering data to pick up the missing 
high-frequency features (sharp edges and jumps) in the recovered conductivity. The method 
applies Total Variation (TV) segmentation and data-driven contrast enhancement to the D- 
bar reconstruction regularized by low-pass filtering with cutoff frequency R. We exploit the 
methodology in mEnu that allows the computation of nonlinear Fourier transform of the 
discontinuous segmented image in the annulus i? — 1 < |A:| < i?, for certain R> R. This new 
transform is added on the annulus R < \k\ < R to the original transform restricted to the 
disc \k\ < i? — 1, and we blend continuously both transforms on the annulus i? — 1 < |A;| < i?. 
From this combined scattering data on the disc |/c| < i?, a new sharper D-bar conductivity 
reconstruction is obtained. This procedure is iterated as outlined in Figure 

Figure [I] shows a reconstruction from simulated FIT data with 0.75% relative noise added 
to the voltage data. Our nonlinear method delivers a piecewise constant and edge-preserving 
reconstruction. The edges are more correctly located near the boundary. This is in accordance 
with the basic intuition about FIT: the deeper in you try to see, the harder it gets. 

Let us comment on the variety of D-bar method we use. The three options are Schrodinger- 
type, 2x2 (first-order) system type and Beltrami-type. They all have two steps: recover 
frequency-domain information and reconstruct via an inverse transform. Theoretically, only 
the Beltrami approach can deal with discontinuous conductivities. However, its second step 
is nonuniform in quality, see [H Section 6.3] and [5l Figure 12]. We use below the shortcut 
method introduced in [5], combining Beltrami-type first step and Schrodinger-type second 
step. 

The rest of this paper is organized as follows. In Section]^ we formulate the mathematical 
FIT model and review the relevant literature. Section [3] is devoted to a discussion of the 
shortcut method. The edge-promoting TV flow is described in Section Section [^introduces 
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Figure 2. Schematic illustration of the nonlinear low-pass filtering approach 
to regularized EIT. The simulated heart-and-lungs phantom (d) gives rise to 
a finite voltage-to-current matrix (orange square), which can be used to 
approximately determine the nonlinear Fourier transform (a). Measurement 
noise causes numerical instabilities in the transform (irregular white patches 
in (a)), leading to an unstable and inaccurate reconstruction (e). However, 
multiplying the transform by the characteristic function of the disc |A;| <5 
yields a lowpass-filtered transform (b), which in turn gives a noise-robust ap¬ 
proximate reconstruction (f). The hybrid method presented in this paper uses 
a priori information about the conductivity to estimate the missing part of the 
nonlinear Fourier transform, resulting in (c). An improved reconstruction (g) 
is achieved. 

a contrast enhancement step based on the CGO sinogram^ which is useful as a robust data- 
fidelity term. The new algorithm we present in this work is tested on two discontinuous 
phantoms for varying levels of noise. Section outlines the computational details and in 
Section the numerical results are presented and discussed. We conclude our findings in 
Section HI 


2. Mathematigal model and literature revievv^ 

We concentrate on the two-dimensional case with Q representing the unit disc. However, all 
our techniques can be extended to simply connected domains 12 C with Lipschitz boundary 
dQ. Throughout the paper the following notation is used. For r >0, 22(0, r) denotes the disc 
in the plane centered at the origin with radius r. For any open set U in the plane, U denotes 
the closure of U. We associate C and by 2 ; = (x, y) — x + iy. 

Let / denote the electric voltage potential maintained at the boundary. The corresponding 
potential u inside the domain 12 satisfies the Dirichlet problem for the elliptic conductivity 
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(a) 


(b) 


(c) 


(d) 


(e) 



Figure 3. Flowchart of the proposed edge-preserving FIT reconstruction method. 


equation 


. X V • a(z)Vu(z) = 0, z e fl 

u\zedii = f{z), z € on, 

where a G L^{^) is an isotropic conductivity satisfying a(z) > c > 0. The Dirichlet-to- 
Neumann (D-N) map, is defined by 


Aaf = 



dQ 


( 2 ) 
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where v denotes the outward facing unit normal vector to the boundary. It is well-known 
that Act : is a bounded linear operator. One can think of A^r as a 

mathematical model for voltage-to-current measurements performed at the boundary. 

The above mathematical model for the inverse conductivity problem was first formulated 
by Alberto Calderon in 1980 m- Inspired by his engineering background, he asked whether 
it is possible to calculate a bounded conductivity a{z) > 5 > 0 G in terms of electrical 

boundary measurements. He gave a proof for the linearized problem in the case of infinite- 
precision knowledge of the noise-free D-N map A^j. 

In practical EIT, one needs to recover the electric conductivity distribution a : ^ M in 

a regularized manner from noisy boundary measurements A^, where ||Acr — A^||y < 5 for a 
known noise level 6 and an appropriate norm H-Hy. Among currently available EIT algorithms, 
the D-bar method is the only one with a regularization analysis m- 

The D-bar method is based on a nonlinear Fourier transform, which is not physically mea¬ 
surable. The definition of the transform depends on certain “almost exponential” functions, 
the so-called complex geometric optics (CGO) solutions. CGO solutions were introduced by 
Faddeev in 1966 [30] and later introduced in the context of inverse problems by Sylvester and 
Uhlmann in 1987 [6^ . 

The D-bar method was invented by Beals and Coifman for use in nonlinear evolution 
equations the early 1980’s, e.g. [9|. The first description of the use of the D-bar approach in 
inverse problems in dimensions n >2 was published by R.G. Novikov in Russian in 1987 (for 
English translation see [58|). The Schrodinger equation approach is still used in dimensions 
three and higher. Rigorous mathematical theory of two-dimensional D-bar methods fall into 
the following three categories. 

Schrodinger-type approach. The first uniqueness result for the inverse conductivity 
problem in dimension two was published by Nachman in 1996 for a G p > 1, in ISZ!- 

The first computational implementation of a D-bar method, based on a kind of Born approxi¬ 
mation, was published in 2000 [63|. The full nonlinear algorithm was developed and equipped 
with a regularization strategy in [saiiEiiiaiizi. For reconstructions from experimental data, 
see jlolESlESlEHlIST]. 

The 2x2 system approach. The smoothness assumption for the conductivity was re¬ 
duced by Brown and Uhlmann to a G p > 2, see m- For further theory and compu¬ 

tational implementation, see usms]. This approach has the benefit of being applicable to 
complex-valued a (conductivity-hi-permittivity), as shown theoretically by Francini [32]. For 
computational results, see [3611351 [37]. 

Beltrami-type approach. Astala and Paivarinta introduced a theoretical reconstruction 
method for a G conductivities in 1311 ]. A corresponding computational reconstruction 
method was developed in OI] and further analysed in |3|. 

Summarizing, numerical D-bar reconstruction methods provide direct, non-iterative, par- 
allelizable, nonlinear reconstruction approaches which are effective on experimental EIT data 
and even in real-time applications [28]. There are also three-dimensional D-bar type methods 

[SHllMlElKEllIIllSe]. 

Conditional stability estimates (involving non-noisy data) for these D-bar methods ap¬ 
peared in mi [7113 E3 ns ED and in dimension greater than two, the most recent stability 
result is in US). For regularization analysis involving noisy data see m- 

Several edge-preserving regularization methods have been suggested for EIT in the liter¬ 
ature, including Total Variation and sparsity-promoting techniques [szmsiisiiEaiMiiis] . 
The method proposed in this paper is very different from all of these approaches as it uses a 
proven regularized D-bar image as its starting point and futher uses the inverse scattering of 
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the D-bar methodology to guide the image segmentation. The proposed method differs from 
[33] . where the CGO sinogram was introduced, by instead focusing on enlarging the scattering 
radius to produce more accurate conductivity reconstructions. 


3. The “shortcut” D-bar reconstruction method 


3.1. The Schrodinger-type D-bar method. The regularized D-bar method [571 l63l ITf] 

for C^-conductivities consists of two steps, namely ^ tR{k) cfr{z). Step 1 is not 
used here; we refer the reader to [53l Chapter 15] for details. 

Step 2 goes from truncated scattering data : C ^ C, supported in the disc |/c| < i?, to 
the regularized conductivity as follows. For each z G fl, solve the integral equation 


(3) 


mR{z,k) = 1 + 


'■"r- i 


tfl(K) 

(27r)2 {k - k)k 


e{—z^ n) mR{z^ n) dRidR2^ 


where e{z^k) := exp {i (fcz -h /cz)} = exp{2i'Re{kz)}. The regularized conductivity is com¬ 
puted by 

(4) aR{z) = {mR{z,0))^ . 


3.2. The Beltrami-type D-bar method. The D-bar method for L^-conductivities laiaii] 
is based on complex geometric optics solutions f±^{z^k) to the Beltrami equation 

(5) dzf±^{z, k) = ±ii{z, k)d^f±^{z, k), 

where /i(z) = (1 — cr(z))/(l -h cf{z)) and f±^{z^k) = e^^^(l -h (9(^)) as | 2 ;| ^ oo. Set 
M±^{z^k) = f±^{z^k). The reconstruction method has three steps: 




Transport 

Matrix 




out of which we only use Step 1 (for Steps 2 and 3, see [53l Section 16.3]). For each fixed 
G C, \k\ < i?, solve 

( 6 ) M±^{;k)\9n + 1= (vi^ + Po)M±^{;k)\9n 

to obtain the CGO traces k) for z G dfl, where and Vo are the projection operators 

described in [1]. 


3.3. The shortcut method. There is a connection [2l Section 5] between the Beltrami 
CGOs of [3111] and Schrodinger CGOs of m, as well as their associated scattering transforms 
r{k) and t(A;), respectively [21 Section 6]. Namely, when cr G we have 

(7) t(A;) = —47Tikr{k) = —2ik f fc) — M-^{z^k))dzidz 2 . 

JM2 

Numerical evidence [5] suggests that the above connection holds for a with jump discontinu¬ 
ities, resulting in numerical equivalence between the reconstructed conductivities. In fact, the 
solution of (© is faster and more stable than the transport matrix method of [1]. Therefore, 
we will use the combined D-bar algorithm presented in [5], called the shortcut method^ which 
has the following steps: 

Step 1 : From noisy boundary measurements to CGO boundary traces by 

solving (§. 
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Step 2: From boundary CGO traces M±^{-,k)\dQ to truncated scattering data tR{k). The 
analyticity of k) outside Q leads to the following development for \z\ > 1 

af(k) a^ik) 

M±^{z, k)^l + + ... 

z z^ 

from which one defines tji{k) := —AmkTR{k), where 


( 8 ) 


TR{k) 


D - «! (fc)) , \k\<R 


\k\ > R. 


Step 3: From truncated scattering data tn^k) to conductivity aR^z). Solve ([s]) and evaluate 


4. Image segmentation method 

In this section, we present a Mumford-Shah based segmentation model. For the simplicity of 
notation, let cro{x) be the input conductivity image defined on 14. The image segmentation 
problem is to find a partition of 14 into K disjoint subdomains {14^1;;}^^, i.e. 


K 

^ ~ U ’ 14;!^ n 14j = 0 for fc ^ j. 

k=i 

The celebrated Mumford-Shah |54| variational segmentation problem is as follows: find 
a piecewise smooth function and a partition edge (closed) set F = d^k such that the 

following functional is minimized: 

-®(/T) = ^ / (fix) - (Toix))‘^dx+ [ \Vf{x)\‘^dx + aR^{T), 

^ Jn Jn\r 

where 44^ (F) denotes the ID Hausdorff measure of the edge set F. Due to the complexity of 
the model, many simplifications are proposed. In the simplest form of the model, / is assumed 
to be piecewise constant on each fhe model is reduced to 


(9) 


mm 

Xk }fc = i 


iC K 

/ {ao{x) - Ckfdx+ '^\dnk\ 

, ^ k=l k=l 


where G M for A: = 1,..., iG is the mean intensity for each subregion The parameter 
A > 0 is used to balance the data fitting and the total length of regions interfaces. This model 
© is hard to be solved directly. In fact, when the regions 14/^ are determined, the optimal Ck 
is given as 


( 10 ) 




Lenu '^ 0 ( 2 ;) 


\^k\ 


Thus we can consider an alternating scheme on solving and iteratively. Once is 
determined, we solve for the By doing so, we introduce the labeling function of the 
disjoint subregions 


( 11 ) 


Uk{x) = 


1 if X ^ Qk 
0 otherwise ’ 


for A: = 1, • • • , K. 
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According to the co-area formula, the perimeter of a set Qk is given by the total variation of 
Uk 

( 12 ) IdQkl = [ \Duk\ 

Jn 

and the total variation \Du\ is defined in distribution sense 

(13) J \Du\ := sup I udiy(^dx : (j) G (7^(fl;M^), |0 (t)| < 1, a.e. x G 

It is well-known that if r G Vk^’^(fl) then \Du\ — \\/u{x)\dx. 

Meanwhile, since each pixel can be only assigned to be one region, the labeling function 
Uk{x) satisfies the following constraint: 

K 

(w) E = 1, a.e. T G 

k=l 


Generally, convex relaxation is made by allowing Ui to take values continuously in [0,1] 
to overcome the computation complexity of the binary constraint. The overall model is 
reformulated as 

. . min f + [ “^^4 

s.t. (ri(t), • • • , uk{x)) G S 


where we denote 
(16) 

the constraint set 


fk{x) = -{(To{x) - Ckf, 


K 


(17) 


S' = < (ri, • • • , uk) ^ BV(fl, [0,1]^) : Uk{x) = 1, a.e. x G 


k=i 


and BV(fl, [0,1]^) denotes the bounded variation functions product space valued on [0,1]^. 

With fixed Ck^ the convex relaxed formulation for Uk allows us to develop efficient algorithms 
based on the well studied total variation minimization. For example, it has been extensively 
studied in [ElEniElEnKIHl. Theoretically, the global solution to the original binary model 
can be achieved for the case of two regions when the intensities Ck are given. 

This above region-based model can be further combined with an edge based approach to 
improve the segmentation quality and speed, such as in [Sana El- Assuming Uk G VF^’^(fl), 
the weighted total variation model can be defined as 

K 

(18) J('u) = V / g{x)\Vuk{x)\dx 

k=iJ^ 


where g{x) >0 is an edge function taking small values at locations with large gradient and 
large values for smooth region. For example, a usual choice is 


(19) 


9{x) 


1 

l + s||Vao(x)||2 


where do is a smoothed version of the given image ao and s > 0 is a positive number. Note 
that if g{x) is identical to 1, it reduces to the model (20). 
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In the following, we present a primal-dual splitting method used in [60l [29l \T7\ [65]. De- 
note u = - ,uk), / = (/i,--- , Ik) and J{u) = Y.k=i Ia 9 {x)\'^Uk\dx and {u, f) = 

fQ^k(x)fk(x)dx, then the convex minimization problem is rewritten as: 

(20) w* = arg min {J('u) + (w,/)} 

ues 

Based on the dual definition ( [T^ , we consider the following min-max model: 

minmax£’(u,p) = {(u,div(p)) -h (u, /)} 
ueS peT 


where p = (pi, • • • ,pk) for pk{x) e pk e {u, div(p)) = Ylk=i In Uk{x)divpkix)dx 

and 

T = {p^ (pi,--- ,Pk) ■■ (^^\Pk{x)\‘^) < g{x), a.e. X en 

1=1 

and for all A; = 1, • • • , |. 


The specific algorithm is given as follows: 


Step 0: Initialization: Choose ci,-- - ,ck as the initial guess of the mean intensity of 
each region. ti,T 2 > 0 are parameters such that tiT 2 <1/8 and G T x S'. 


^0 = ^0 


Step 1: 
Step 

Step 

Step 

Step 


Set i := 0 and run the outer loop as follows: 

Set j := 0 and run the inner loop to eompute Uk{x) for A: = 1, • • • and t G 
1.1: Compute the dual variable — By (p-^ +ri'Vu^) where By (•) denotes the 
projection operator onto the convex set T. 

1 . 2 : Compute the primal variable = B^(i6-^ +T 2 (/^ + div*jy<^)) where B^(-) 
denotes the projection operator onto the convex set S. 

1.3: Compute the auxiliary primal variable: — u^) 


1.4: Set j — j + 1 update until stopping conditions satisfied, output u. 

Step 2: Compute the piecewise regions by the binarization of Uk{x). Generally, a global 


( 21 ) 


minimizer of (20) might not be binary, and a final thresholding step needs to be 
taken to get a binary solution 

X _ f 1 ul{x) =v[idQ^{ui{x),U2{x),... ,uk{x)} 

^k\^) I Q otherwise 

If the maximizer is not unique, the maximizer with smallest subscript is used as 
a convention. 


Update the mean intensity estimation: c\ for k — 1, • • • by (10) and by 

(P. 


Step 3: 

Step 4: Set i = z + 1 update until stopping conditions satisfied. 


After we obtain the label functions and for A; = 1, • • • , K, an image can be 

reconstructed as a piecewise constant function with the mean intensity in the corresponding 
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k— th region, i.e. 
( 22 ) 


K 

o-Tv(x) = '^ul(x)ck 
k=l 


Hence, this reconstructed image can be used as a piecewise constant regularized approximation 
to the original image by letting ao(x) := (Tdb{x) in the whole algorithm described in Figure 

El 


5. Data-Driven Contrast Adjustments 

5.1. The Beltrami CGO sinogram. We extend the concept of the CGO sinogram, intro¬ 
duced in [22]) to discontinuous conductivities. Set 

Sa{0,^,p) ■■= - 1 , 

where 2 : = and k = for ^ [—7r,7r) and the traces of the CGO 

solutions are solved from the noisy FIT data using equation The radius p must be 
smaller than the noise-dependent cutoff frequency R. 

The traditional data-fidelity term used in EIT is ||Acr — Acr/||. We use instead the CGO 
sinogram data-fidelity term 

(23) ||5cr(0, p) — p) II ]^2(^j2y 

where denotes the two-dimensional torus. 


5.2. Contrast enhancement. Assume that the piecewise constant conductivity satisfies 
supp((j — 1) C ri and that we know a priori approximate bounds 0 < c < 1 and C > 1 
such that 


(24) 


min a(z) > c, msixa(z) < C. 

zen zeQ 


Let a denote an approximate reconstruction to the true conductivity a defined on 12, whose 
contrast we intend to improve, and suppose supp(a — 1) C 12. Set f{z)=a — l and denote 


(25) 


m = min/( 2 ;), 

Z^il 


M = max f{z), 

Z^il 


and assume that m < 0 and M > 0. Let s and t be two parameters such that 0 < 5 < 1 and 
0 < t < 1 and define 

( t{C — l)f{z)/M for z satisfying n{z) > 1, 

(26) := 1 + < 5(c — l)/( 2 ;)/m for z satisfying n{z) < 1, 

[ 0 otherwise. 

Note that (Jo,o = 1? cind the maximum values 5 = 1 and t = 1 yield the maximal image 
contrast below and above 1, respectively, i.e. 


mincTi Az) = c, maxcr. Az) = C. 


We determine the optimal values for as the minimizers of the nonlinear data-discrepancy 
functional based on the CGO sinogram: 


(27) 


(so, to) arg min 

(sj)e[o,i]2 


(•, • )P) ~ , • )P)IIl2(t2) 

, • ,p)||L2(Tr2) 


The result of contrast-enhancement is then acE •— ^so2o- 
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cri 

C^2 

Background 1.0 
Lung 0.5 

Heart 2.0 

Pipe 1.0 

Top layer (oil) 1.2 

Middle layer (water) 2.0 
Bottom layer (sand) 0.3 


Table 1. Conductivity values in the two simulated phantoms shown in Figure 
Left: heart-and-lungs phantom ai. Right: oil pipeline phantom (J 2 . 



Figure 4. True conductivity phantoms; for conductivity values see Table 
Left: heart-and-lungs phantom ai. Right: pipeline phantom (72. 


In this paper, the objective function in (27) is minimized via the DIRECT algorithm 
in an analogous fashion to 


The expression DIRECT refers to “Dividing RECTangles” 
which suggests the strategy of this sampling, global search algorithm. 


6. Numerigal Implementation 

6.1. Simulation of noisy EIT data. Our numerical experiments deal with two simulated 
discontinuous conductivity phantoms, namely a heart-and-lungs phantom ai and a cross- 
section of a stratified oil pipeline phantom (J 2 . See Tableand Figure]^ 

6.1.1. Computation of the discrete Neumann-to-Dirichlet map. EIT data was simulated us¬ 
ing the finite element method and following [53l Sections 13.2.3 and 16.3.3]. We use the 
trigonometric basis functions 

... f cos((n-h l)6^/2) for odd n, 

9n[d)^\ 1/2^ 


TT ^/^sin(n0/2) 


for even n. 


where 1 < n < 2N. The Neumann-to-Dirichlet map TZ^ is approximated by the matrix 
Ro- = [(Ra)m,n] given by 


rZTT 

(r^cr)m,n ~ (1^(7 0715 0 m) ~ (^) 0 m(^) ^^5 

Jo 

where Rcr0n = Un\dQ with V •a'Vun = 0 in D, (a {dun/dv'))\dQ = 0^ and dS — 0. Here, 

1 < m,n < 2N. In this paper, we use = 16 which corresponds to 33 linearly independent 
current patterns. 
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Relative Gaussian noise was added to the boundary voltage data as in [M]. Namely, define 
with (Rcr) 77 T,,n — where 

(28) ~ “ 1 “ V || || j^oo , 

T] denotes the noise level so that 10077% noise is added, and A/i,...,A/ 2 W+i are independent 
Gaussian distributions with mean zero and variance one, which are implemented through the 
MATLAB function randn generating pseudorandom values drawn from the standard normal 
distribution A/’(0,1). The noisy D-N matrix is then formed by inverting Rcr and adding a zero 
row and column on top and left as described in [53] . 

The algorithm outlined in Figure]^ was applied to the two phantoms shown in Figure]^ 
In this paper, the algorithm was performed for J = 3 iterations, and for three noise levels as 
follows: zero added noise, 0.1% added noise and 0.75% added noise (corresponding to setting 
77 = 0, 0.001, 0.0075). The output of the algorithm after the J = 3 iterations is denoted by 

^CE • 

6 . 2 . Computational grids. All the computations on the 2 ;-plane (D-bar reconstructions, 
TV flow outputs and CE outputs) were generated on a z-grid of 2 ^x 2 ^ equidistributed points 
of the square [—s, s) x [—s, s) with £ = 8 and s = 2.3. Thus, the z-grid consists of 2^^ = 65536 
points. The D-bar equation solver (used to solve Q) was executed just on the set of 9729 
points of the z-grid belonging to the closed disc Q. Note that the larger z-region is needed to 
extend the scattering data. 

The k-ghds are problem specific, i.e. for lower levels of noise a larger radius i? can be 
used for the initial low-pass filtering in the nonlinear Fourier domain. In each case, we fixed 
the two parameters i? (the initial low pass filtering chosen intuitively by looking at where 
the scattering data “blows up” in magnitude), and i? (the increased scattering radius to be 
determined by solving the Beltrami equation). In each case, the fc-grid for the scattering data 
was comprised of 2^ x 2^ equispaced points on the square [—i?, i?) x [—i?, i?). The traditional 
(and original) D-bar image is computed from scattering data satisfying \k\ < R and all 
subsequent D-bar images computed using the larger disc |A:| < i?. 

Remark: If supp((j — (Jq) C D(0, 1) for some positive constant (Jq % 1, the algorithm can be 
re-scaled as follows. Defining a := cr/ao, we have supp(a — 1) C D(0,1). Apply the above 

algorithm to = croAcr and write for the output generated. Take (Jq o^ce* final 

approximation to a. 


6.3. Computation of the initial nonlinear Fourier transform. This step corresponds 
to Figure l^a). From the boundary measurements A^ we use and Q to compute the 
initial scattering data r^{k) on a fc-disc of radius R. We write t^(A:) := —47vikr^{k). 

6.4. D-bar reconst ruc tion. This step corresponds to Figure [^b). Use the shortcut method 

to compute from the scattering transform U“^(fc). If j = 1 then 


3.3 


described in Section 
use the smaller cutoff disc \k\ < R^ else use larger disc |A:| < i?. 


6 . 


5. Edge-enhancement using TV flow. This step corresponds to FigureJ^c). Introduce 

- - - - - - Aj) - - - - . „ „ „ - rV_. 


edges into the D-bar image cr^g by applying the segmentation flow of Section^ The resulting 


piecewise constant image, defined by (|22|), is called Note that the initial guess of the mean 


intensity are obtained directly by K-means algorithm, where K denote the number regions 
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pre-selected. In practice, this is a reasonable guess at the number of regions of different 
conductivity in your domain. 


6.6. Contrast enhancement. This step corresponds to Figure [^d). The CGO sinograms 
were implemented via 33x33 matrices [(S)^^/] with ki) — l = 26^^^^) — 

1, where 

(9^ = (m - 1 - N)2ti/{2N + 1), l<m<2N + 1, 


with N — 16, — 01^ and M^{z^k) refers to the solution explained in Section 5.1 for both 

the true a and the corresponding approximations Therefore, {9m} and {ipi} are the same 
partition of the interval (—tt, tt). Note that one can choose the points ipi independently of the 
9 values if desired. 

Finally, the output image is determined by plugging a — into (26) for (so,to) G 
[0,1]^ obtained via the DIRECT optimization strategy (see 


6.7. Extension of the Scattering Transform. This step corresponds to Figure [^e). The 
radius of the admissible scattering data is increased from i? to i? by computing new stable 
scattering data in the annulus i? — 1 < |fc| < i? corresponding to the TV-sharpened and 
contrast adjusted image follows. First, evaluate the Beltrami coefficient fi{z) = (1 — 

^ce(^))/(1+ ^ce(^)) on the z-grid [—2.3, 2.3) x [—2.3, 2.3). Next, solve the Beltrami equation 
^ for the CGO solutions f±ii{z^ k) for k in fc-annulus i? — 1 < |A;| < i?. Finally, evaluate the 
scattering data T{k) in the fc-annulus i? — 1 < |A:| < i? via 

(29) G /■ - M.^{z,k)) dzidz2 

where M±ij^{z^k) = f±^{z^k). Call this new scattering data f^^\k). The new scattering 

data on the larger radius is then 

T^^\k) ■= x(fc)r°(/c) + (1 - x{k)) A)(fc), |fc| < R, 

where we used the polynomial radial cutoff function x defined by 

[ 1, if |A:| <R-1, 

= l p{\k\-{R-l)), iiR-l<\k\<R, 

[ 0, if |A;| > i?, 

with p{t) = 1 — 3t^ -h 2t^ to blend the data in the overlap region i? — 1 < |A:| < i?. Note that 

in Figure]^ we use the notation P{k) := —AnikT^^^k) and (k) := —4:7Tikf^^\k). 


6.7.1. Computation of Scattering Data. We need a comparison for our new extended 

scattering data with the best possible scattering data. By best possible scattering data we 
refer to the scattering data that is obtained by computing the CGO solutions to the Beltrami 
equation i) with fi corresponding to the true cr, and evaluating the scattering data r{k) via 
(29). For details on how to solve the Beltrami equation and generate the scattering data r, 
the reader is referred to mum- 


7. Numerical Results 

The algorithm was tested on the two phantoms shown in Figure and the results are 
shown here. 
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7.1. Example 1: A Heart and Lungs Phantom. For the heart and lungs phantom, (Ji, 
it was assumed known apriori that the internal conductivity was bounded between c = 0.3 
and c = 2.5. Setting the initial scattering radius R was 5 and the enlarged radius R to 10 
for all noise levels proved sufficient. The initial scattering data was reliable for all noise levels 
within the fc-disc of radius 5. The parameters for the TV flow were K = A and A = 0.1. 

The scattering data for each noise level is displayed in Figure The figures contain images 
of the actual Beltrami scattering transform tb on the larger k disc of radius 10. This is used 
to evaluate the efficacy of the new proposed approach. The true scattering data tb was 
computed by solving (§ and ( [2^ with the known /i = serve as a best ease 

seenario baseline. 


Re 

Im 

Re 

Im 

Re 

Im 

% 

"■b * y ' 


« 
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^ ^ • 


a) no added 

[ noise 

b) 0.1% 

noise 

c) 0.75% 

noise 



- 0.15 - 0.1 - 0.05 0 0.05 0.1 0.15 


Figure 5. Images of the real and imaginary parts of the reliable scattering 
transform tb (computed directly from the Beltrami equation) and the com¬ 
bined scattering data for the two iterations and from the simulated 
Dirichlet-to-Neumann corresponding to the heart and lungs phantom ai. 

The reconstructed conductivities for each stage of the algorithm are displayed in Figure]^ 
Note that the reconstructions are displayed on the same color scale as the original conductivity 
shown in Figure]^ (left) for ease of comparison. 

Tables [^1^ and show the relative errors and the Structural SiMilarity Index (SSIM) 
values for the heart-and-lungs phantom ai of Example 1 for zero added noise, 0.1% added 
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True (Ji 

<0 


Iter 1 

Iter 2 

Iter 3 


Figure 6. This figure shows the real parts of the numerical approximations 
obtained for the heart and lungs phantom of Example 1, i.e. ai. The pic¬ 
ture consists of three parts a), b), c) corresponding to our three cases of added 
simulated noise in the EIT voltage data: zero added noise, noise of relative am¬ 
plitude 0.1% and noise of relative amplitude 0.75%, respectively. For the sake 
of comparison, the true conductivity ai is displayed above the reconstructions 
(all on the same color scale). For each noise level, the D-bar reconstruction 
^DB columns), the TV sharpened image (middle columns), and the 

(j\ 

contrast adjusted TV sharpened images The first, second, and third rows 
correspond to the first, second and third iterations, respectively. 


^DB CTxv CTCE (Jtv ^CE 


C^DB ^TV ^CE 


i* i <i§ t*f # (I*! /I 

i§^§^§ ^§'i§^§ €%«>•> 


a) no added noise 


b) 0.1% noise 


c) 0.75% noise 


noise, and 0.75% added noise, respectively. The error values are presented for each step of the 

ii) ii) 

proposed algorithm: the D-bar reconstruction the sharpened reconstruction cr.fy, and 

(j\ 

the contrast adjusted sharpened reconstruction 

Relative Error SSIM 


3 

"dr 

CTtv 

^CE 

3 

"nR 

CTtv 

^CE 

1 

0.1240 

0.1240 

0.1202 

1 

0.6600 

0.6600 

0.6541 

2 

0.1095 

0.1185 

0.1168 

2 

0.7351 

0.6903 

0.6878 

3 

0.1054 

0.1157 

0.1145 

3 

0.7425 

0.7117 

0.7096 


Table 2. The relative errors as well as the SSIM values for the zero added 
noise case for Example 1 with the heart and lungs phantom ai. 


In the case of zero and 0.1% added relative noise, the relative error decreases with each 
iteration and the SSIM value increases. Both measures confirm that the image is “improving”. 








16 


S. J. HAMILTON, J. M. REYES, S. SILTANEN, AND X. ZHANG 


Relative Error 


SSIM 


Table 3. 


j 


CTtv 

'^CF, 

j 


CTtv 

'^CF. 

1 

0.1009 

0.1233 

0.1194 

1 

0.7304 

0.6603 

0.6545 

2 

0.1083 

0.1174 

0.1158 

2 

0.7348 

0.6907 

0.6883 

3 

0.1009 

0.1142 

0.1133 

3 

0.7304 

0.7131 

0.7111 


The relative errors as well as t 


le SSIM values for the 


0 .1% added 


noise case for Example 1 with the heart and lungs phantom ai. 


Relative Error SSIM 


j 

‘^RR 

cr^V 

'^CF 

j 

‘^RR 

cr^V 

'^CF 

1 

0.1092 

0.1202 

0.1167 

1 

0.6897 

0.6680 

0.6639 

2 

0.1076 

0.1165 

0.1134 

2 

0.6964 

0.6816 

0.6776 

3 

0.1092 

0.1151 

0.1171 

3 

0.6897 

0.6773 

0.6790 


Table 4. The L? relative errors as well as the SSIM values for the 0.75% 
added noise case for Example 1 with the heart and lungs phantom ai. 


In the case of 0.75% added relative noise, the L? relative error and SSIM remain approximately 
the same but visually one can see the improvements in each iteration as the artifact present 
in iterations 1 and 2 (above the heart) is absent in iteration 3. 

7.2. Example 2: An Industrial Pipeline. For the cross-section of the industrial pipeline 
phantom, (J 2 , it was assumed known apriori that the internal conductivity was bounded 
between c = 0.1 and c = 2.5. Here we used The parameters for the TV flow were A = 5 in 
the TV flow and Table gives the values of the scattering radii and the A parameter used in 
for the TV flow for each noise level. 


Added Noise Level 

0 % 

0 .1% 

0.75% 

R 

6 

5 

4 

R 

10 

8.3 

6.6 

X 

0.3 

0.5 

0.5 


Table 5. Parameter values for the new algorithm for each noise level for the 
industrial pipeline phantom (J 2 . 


The scattering data for each noise level is displayed in Figure Here the figures con¬ 
tain images of the actual Beltrami scattering transform tb on the larger k disc of radii of 
10, 8.3, 6.6 for the varying noise levels respectively. This is used to evaluate the efficacy of 
the new proposed approach. The true scattering data tb was computed by solving and 
(29) with the known /i = serve as a best ease seenario baseline. 

The reconstructed conductivities for each stage of the algorithm are displayed in Figure]^ 
Note that the reconstructions are displayed on the same color scale as the original conductivity 
shown in Figure]^ (right) for ease of comparison. Furthermore, note that the ring of constant 
conductivity along the boundary (representing the thickness of the pipe) has been enforced 
in the reconstructions as this can also be considered apriori information for this application. 

Tables and show the relative errors and the Structural SiMilarity Index (SSIM) 
values for the industrial pipe phantom a 2 of Example 2 for zero added noise, 0.1% added 
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Figure 7. Images of the real and imaginary parts of the reliable scattering 
transform tb (computed directly from the Beltrami equation) and the com¬ 
bined scattering data for the two iterations and from the simulated 
Dirichlet-to-Neumann corresponding to the industrial pipeline phantom (J 2 . 

Left: zero added noise and ^ = 10, Middle: 0.1% added noise and ^ = 8.3, 

Right: 0.75% added noise and R = 6.6. 

noise, and 0.75% added noise, respectively. The error values are presented for each step of the 

proposed algorithm: the D-bar reconstruction the sharpened reconstruction cr^fy, and 

(i) 

the contrast adjusted sharpened reconstruction 

Relative Error SSIM 


3 


^TV 

'^CF, 

j 

‘^nR 

^TV 

'^CF. 

1 

0.1926 

0.1926 

0.2157 

1 

0.7097 

0.7097 

0.7130 

2 

0.1813 

0.1985 

0.2286 

2 

0.7127 

0.6976 

0.6972 

3 

0.1862 

0.2069 

0.2228 

3 

0.7054 

0.6130 

0.6128 


Table 6 . The relative errors as well as the SSIM values for the zero added 
noise case for Example 2 with the pipeline phantom (J 2 . 


In the zero added noise case we see that the relative error in the D-bar images decreases 
with each iteration and the SSIM stays approximately the same. Interestingly, the sharpened 
and contrast adjusted images appear to perform slightly worse under the metrics. However, 
visually the sharpened images much better reflect the physical scenario (oil, water, sand) then 
their smooth D-bar counterparts as they contain nice clean divisions between the layers. As 
the noise level increases, the method is still able to clearly distinguish between oil, water and 
sand at even the first iteration and thus the algorithm could be stopped there (i.e. at cr^y) 
for each noise level. 
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True (72 



CTxv CTCE (Jtv ^CE ^TV ^CE 



a) no added noise b) 0.1% noise 



i 


Figure 8 . This figure shows the real parts of the numerical approximations 
obtained for the industrial pipe phantom of Example 2, i.e. (J 2 . The pic¬ 
ture consists of three parts a), b), c) corresponding to our three cases of added 
simulated noise in the FIT voltage data: zero added noise, noise of relative am¬ 
plitude 0.1% and noise of relative amplitude 0.75%, respectively. For the sake 
of comparison, the true conductivity (72 is displayed above the reconstructions 
(all on the same color scale). For each noise level, the D-bar reconstruction 

(Jbb (l^ft columns), the TV sharpened image a^y (middle columns), and the 

(i) 

contrast adjusted TV sharpened images The first, second, and third rows 
correspond to the first, second and third iterations, respectively. 


Relative Error 


SSIM 


j 


CTtv 

U) 

^CE 

j 

‘^DB 

CTtv 

^CR 

1 

0.1970 

0.2176 

0.2479 

1 

0.7292 

0.6834 

0.6829 

2 

0.2279 

0.2309 

0.2305 

2 

0.6982 

0.6826 

0.6837 

3 

0.2329 

0.2404 

0.2392 

3 

0.6997 

0.5744 

0.5653 


Table 7. The relative errors as well as t 


le SSIM values for the 0.1% added 


noise case for Example 2 with the pipeline phantom (J 2 . 


Regarding the contrast enhanced images, recall that the approximate upper and lower 
bounds C = 2.5 and c = 0.1, respectively, were used. In practice, closer approximations 
may be known (in particular for the case of oil, water, and sand) therefore improving the 
reconstructed values. An investigation into the optimal parameters c and C and minimization 
scheme for the contrast adjustment is outside the scope of this introductory paper. 
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Relative Error SSIM 


j 


CTtv 

'^CF, 

j 


CTtv 

'^CF. 

1 

0.2043 

0.2181 

0.2617 

1 

0.7123 

0.6908 

0.6958 

2 

0.2255 

0.2273 

0.2434 

2 

0.6873 

0.6972 

0.6960 

3 

0.2604 

0.2758 

0.2622 

3 

0.6933 

0.5868 

0.5690 


Table 8 . The L? relative errors as well as the SSIM values for the 0.75% 
added noise case for Example 2 with the pipeline phantom (72. 


8. Conclusions 

EIT data contains information about the conductivity in an indirect, nonlinear and unsta¬ 
ble way. Theoretically, [3] shows that infinite-precision data (the DN map) contains enough 
information to uniquely determine the conductivity. However, a practical data matrix is a 
noisy and finite-dimensional approximation of infinite-dimensional data A^r, and in most cases 
does not actually correspond to any conductivity (as A^ is not in the range of the forward map 
(7 i-A Act). Therefore, EIT reconstruction methods need to be regularized to yield noise-robust 
results. Regularization is based on complementing the insufficient measurement data by a 
priori information about the conductivity. 

Currently there are not many regularized reconstruction methods for EIT. The theory of 
Tikhonov regularization and related variational methods applies to a wide class of nonlinear 
forward maps [Hj, but, alas, not to the extremely nonlinear case of EIT. For partial results, 
see ISOl sa HD. The enclosure method for detecting convex hulls of inclusions admits a 
regularization analysis [39], but it only yields partial information (e.g. information about the 
locations of inclusions rather than their conductivity values). The D-bar method [47] is the 
only regularized reconstruction method that produces actual conductivity images, but the 
reconstructions are always smooth because of a nonlinear low-pass filter involved. Indeed, 
the assumptions of the regularized D-bar method include continuous differentiability of the 
conductivity and thus the smoothing is not unexpected. 

In many applications of EIT, such as nondestructive testing, the conductivity distribution 
can be assumed to be piecewise constant. This is approximately the case in medical imaging 
as well. Therefore, it is desirable to design a regularized reconstruction method producing 
piecewise constant images. 

In this paper, a noise-robust EIT reconstruction method that always results in a piecewise 
constant image was both presented and tested on simulated noisy EIT data. Therefore, this 
paper demonstrates that one can achieve the above goal, at least partially. The authors 
note that this is an initial feasibility study only, and do not prove that the new combined 
method is itself a regularization strategy. However, there may be hope to prove that the 
reconstruction approaches the true piece-wise constant conductivity along a stable path as 
the data error tends to zero. Namely, the low frequencies of the reconstruction are provided 
by the regularized D-bar method, while the large frequencies are built based on the piecewise 
constant assumption. The limit frequencies between the different treatments can be assumed 
to tend to infinity as the noise level vanishes. 

Although the scattering data in the extended annuli for and in each example is not 
identical to that of tb, nonetheless, the subsequent corresponding conductivity reconstructions 
show marked improvements in the locations, sharpness of edges, and conductivity values of 
the inclusions. Take particular note of the 0.75% noise case of Example 1 where in iterations 
1 and 2 the images (7 db, ^ltv, and acE all contain a strong artifact above the heart which 
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is not present after an additional iteration (i.e. in iteration 3). In Example 2, we see that 
after a single sharpening of the original D-bar image the corresponding is enough for the 
goal of distinguishing between oil, water, and sand in the pipeline. These two cases (shown 
for various levels of noise) demonstrate the flexibility of the algorithm: for some cases it is 
appropriate to be applied iteratively for improvements and to remove artifacts, whereas in 
other cases a single sharpening iteration is adequate. 

Additional modifications to the proposed method can be easily applied. In particular, in 
lieu of fixing the maximal number of iterations J beforehand, alternative stopping criteria 
could be applied. E.g., for j > 1 instead of asking above if j — J, return fhe final 

image if 


Jj) Jj-i) 



^CE '^CE 

P ^ 



P 


< thresh 


for some specified threshold thresh. Another option concerns the radii R and R. In the 
examples presented here, the radii R and R are fixed throughout the algorithm for each data 
set Act. An alternative approach could compute the Beltrami scattering data on progressively 
larger annuli so that in Step (e) of the algorithm in Figure]^ the new scattering data r^^\k) 
is computed for i? — 1 < |A:| < R + {j — l)Ai? given a fixed stepsize AR > 0 (j > 1). Such 
approaches, while interesting are outside the scope of this work. 


Acknowledgments 

S. J. Hamilton, J. M. Reyes, and S. Siltanen were supported by the Academy of Finland 
(Finnish Centre of Excellence in Inverse Problems Research 20122017, decision number 250215) 
S. J. Hamilton was additionally supported by Sal We Research Program for Mind and Body 
(Tekes - the Finnish Funding Agency for Technology and Innovation grant 1104/10). J. M. 
Reyes was additionally supported by the Engineering and Physical Sciences Research Coun¬ 
cil (EPSRC), reference EP/K024078/1. X. Zhang was supported by NSFC11101277 and 
NSFC91330102. 


References 

[1] K. Astala, J.L. Mueller, L. Paivarinta, A. Peramaki, and S. Siltanen. Direct electrical impedance tomog¬ 
raphy for nonsmooth conductivities. Inverse Problems and Imaging, 5(3):531-549, 2011. 

[2] K. Astala, J.L. Mueller, L. Paivarinta, and S. Siltanen. Numerical computation of complex geometrical 
optics solutions to the conductivity equation. Applied and Computational Harmonie Analysis, 29(1):391- 
403, 2010. 

[3] K. Astala and L. Paivarinta. A boundary integral equation for Calderon’s inverse conductivity problem. 
In Proe. 7th Internat. Conferenee on Harmonie Analysis, Colleetanea Mathematiea, 2006. 

[4] K. Astala and L. Paivarinta. Calderon’s inverse conductivity problem in the plane. Annals of Mathematies, 
163(l):265-299, 2006. 

[5] K. Astala, L. Paivarinta, J. M. Reyes, and S. Siltanen. Nonlinear Fourier analysis for discontinuous 
conductivities: Computational results. Journal of Computational Physies, 276:74-91, 2014. 

[6] E. Bae, J. Yuan, and X.C. Tai. Global minimization for continuous multiphase partitioning problems using 
a dual approach. International journal of eomputer vision, 92(1):112-129, 2011. 

[7] J. A. Barcelo, T. Barcelo, and A. Ruiz. Stability of the inverse conductivity problem in the plane for less 
regular conductivities. Journal of Differential Equations, 173(2):231-270, 2001. 

[8] T. Barcelo, D. Faraco, and A. Ruiz. Stability of Calderon inverse conductivity problem in the plane. 
Journal de Mathematiques Pures et Appliques, 88(6):522-556, 2007. 

[9] Richard Beals and Ronald R. Coifman. Multidimensional inverse scatterings and nonlinear partial differ¬ 
ential equations. In Pseudodifferential operators and applieations (Notre Dame, Ind., 1984), pages 45-70. 
Amer. Math. Soc., Providence, RI, 1985. 










HYBRID SEGMENTATION AND D-BAR METHOD FOR EIT 


21 


[10] E. Beretta and E. Erancini. Lipschitz stability for the electrical impedance tomography problem: the 
complex case. ArXiv e-prints^ August 2010. 

[11] J. Bikowski, K. Knudsen, and J. L. Mueller. Direct numerical reconstruction of conductivities in three 
dimensions using scattering transforms. Inverse Problems, 27:19pp, 2011. 

[12] G. Boverman, D. Isaacson, T.-J. Kao, Saulnier, G. J., and J. G. Newell. Methods for direct image recon¬ 
struction for EIT in two and three dimensions. In Proeeedings of the 2008 Eleetrieal Impedanee Tomography 
Conferenee, Dartmouth Gollege, in Hanover, New Hampshire, USA, June 16 to 18 2008. 

[13] X. Bresson, S. Esedoglu, P. Vandergheynst, J.P. Thiran, and S. Osher. East global minimization of the 
active contour/snake model. Journal of Mathematieal Imaging and Vision, 28(2):151-167, 2007. 

[14] R. M. Brown and G. Uhlmann. Uniqueness in the inverse conductivity problem for nonsmooth conductiv¬ 
ities in two dimensions. Communieations in Partial Differential Equations, 22(5): 1009-1027, 1997. 

[15] A.-P. Galderon. On an inverse boundary value problem. In Seminar on Numerieal Analysis and its Ap- 
plieations to Continuum Physies (Rio de Janeiro, 1980), pages 65-73. Soc. Brasil. Mat., Rio de Janeiro, 
1980. 

[16] P. Garo, A. Garcia, and J. M. Reyes. Stability of the Galderon problem for less regular conductivities. 
Journal of Differential Equations, 254:469-492, 2013. 

[17] A. Ghambolle and T. Pock. A first-order primal-dual algorithm for convex problems with applications to 
imaging. Journal of Mathematieal Imaging and Vision, 40(1):120-145, 2011. 

[18] Antonin Ghambolle, Stacey E. Levine, and Bradley J. Lucier. An upwind hnite-difference method for total 
variation-based image smoothing. SIAM J. Imaging Set, 4(l):277-299, 2011. 

[19] Tony E Ghan and Selim Esedoglu. Aspects of total variation regularized 1 1 function approximation. SIAM 
Journal on Applied Mathematies, 65(5): 1817-1837, 2005. 

[20] Tony E. Ghan, Selim Esedoglu, and Mila Nikolova. Algorithms for finding global minimizers of image 
segmentation and denoising models. SIAM Journal of Applied Mathematies, 66(5):1632-1648, 2006. 

[21] M. Gheney, D. Isaacson, and J. G. Newell. Electrical impedance tomography. SIAM Review, 41(1):85-101, 
1999. 

[22] Eric T. Ghung, Tony E. Ghan, and Xue-Gheng Tai. Electrical impedance tomography using level set 
representation and total variational regularization. Journal of Computational Physies, 205(1):357 - 372, 
2005. 

[23] A. Glop, D. Earaco, and A. Ruiz. Stability of Galderon’s inverse conductivity problem in the plane for 
discontinuous conductivites. Inverse Problems and Imaging, 4(1):49-91, 2010. 

[24] H. Gornean and K. Knudsen. Reconstruction from one boundary measurement of a potential homoge¬ 
neous of degree zero. Journal of Inverse and Ill-Posed Problems, 13(3-6) :413-423, 2005. Inverse problems: 
modeling and simulation. 

[25] M. DeAngelo and J. L. Mueller. 2d D-bar reconstructions of human chest and tank data using an improved 
approximation to the scattering transform. Physiologieal Measurement, 31:221-232, 2010. 

[26] Eabrice Delbary and Kim Knudsen. Numerical nonlinear complex geometrical optics algorithm for the 3D 
Galderon problem. Inverse Problems and Imaging, 8(4):991-1012, 2014. 

[27] D. G. Dobson and E. Santosa. An image enhancement technique for electrical impedance tomography. 
Inverse Problems, 10:317-334, 1994. 

[28] Melody Dodd and Jennifer L. Mueller. A real-time d-bar algorithm for 2-d electrical impedance tomography 
data. Inverse Problems and Imaging, 8(4): 1013-1031, 2014. 

[29] Tony E. Ghan Ernie Esser, Xiaoqun Zhang. A general framework for a class of first order primal-dual 
algorithms for convex optimization in imaging science. SIAM J. Imaging Set, 3:1015-1046, 2010. 

[30] L. D. Eaddeev. Increasing solutions of the Schrodinger equation. Soviet Physies Doklady, 10:1033-1035, 
1966. 

[31] D. Earaco and K. Rogers. The Sobolev norm of characteristic functions with applications to the Galderon 
Inverse Problem. The Quarterly Journal of Mathematies, 64:133-147, 2013. 

[32] E. Erancini. Recovering a complex coefficient in a planar domain from Dirichlet-to-Neumann map. Inverse 
Problems, 16:107-119, 2000. 

[33] S. J. Hamilton, A. Hauptmann, and S. Siltanen. A Data-Driven Edge-Preserving D-bar Method for Elec¬ 
trical Impedance Tomography. Inverse Problems and Imaging, 8(4): 1053-1072, 2014. 

[34] S. J. Hamilton, M. Lassas, and S. Siltanen. A Direct Reconstruction Method for Anisotropic Electrical 
Impedance Tomography. Inverse Problems, 30:(075007), 2014. 

[35] Sarah J Hamilton and Jennifer L Mueller. Direct eit reconstructions of complex admittivities on a chest¬ 
shaped domain in 2-d. IEEE transaetions on medieal imaging, 32(4):757-769, 2013. 



22 


S. J. HAMILTON, J. M. REYES, S. SILTANEN, AND X. ZHANG 


[36] S.J. Hamilton, C.N.L. Herrera, J. L. Mueller, and A. VonHerrmann. A direct D-bar reconstruction algo¬ 
rithm for recovering a complex conductivity in 2-D. Inverse Problems, 28:(095005), 2012. 

[37] C.N.L. Herrera, M.F.M. Vallejo, J.L. Mueller, and R.G. Lima. Direct 2-d reconstructions of conductivity 
and permittivity from eit data on a human chest. Medieal Imaging, IEEE Transaetions on, 34(l):267-274, 
Jan 2015. 

[38] M. Huhtanen and A. Peramaki. Numerical solution of the R-linear Beltrami equation. Mathematies of 
Computation, 81:387-397, 2012. 

[39] M. Ikehata and S. Siltanen. Electrical impedance tomography and Mittag-Leffler’s function. Inverse Prob¬ 
lems, 20:1325-1348, 2004. 

[40] D. Isaacson, J.L. Mueller, J.C. Newell, and S. Siltanen. Imaging cardiac activity by the D-bar method for 
electrical impedance tomography. Physiologieal Measurement, 27:S43-S50, 2006. 

[41] Bangti Jin and Peter Maass. An analysis of electrical impedance tomography with applications to Tikhonov 
regularization. ESAIM Control Optim. Cale. Var., 18(4):1027-1048, 2012. 

[42] Bangti Jin and Peter Maass. Sparsity regularization for parameter identification problems. Inverse Prob¬ 
lems, 28(12):123001, 70, 2012. 

[43] J.P. Kaipio, V. Kolehmainen, E. Somersal, and M. Vauhkonen. Statistical inversion and monte carlo 
sampling methods in electrical impedance tomography. Inverse Problems, 16(5): 1487-1522, 2000. 

[44] B. Kaltenbacher, A. Neubauer, and O. Scherzer. Iterative regularization methods for nonlinear ill-posed 
problems, volume 6. de Gruyter, 2008. 

[45] K. Knudsen. On the Inverse Conduetivity Problem. PhD thesis. Department of Mathematical Sciences, 
Aalborg University, Denmark, 2002. 

[46] K. Knudsen, M. Lassas, J.L. Mueller, and S. Siltanen. D-bar method for electrical impedance tomography 
with discontinuous conductivities. SIAM Journal on Applied Mathematies, 67(3) :893, 2007. 

[47] K. Knudsen, M. Lassas, J.L. Mueller, and S. Siltanen. Regularized D-bar method for the inverse conduc¬ 
tivity problem. Inverse Problems and Imaging, 3(4):599-624, 2009. 

[48] K. Knudsen, J.L. Mueller, and S. Siltanen. Numerical solution method for the dbar-equation in the plane. 
Journal of Computational Physies, 198:500-517, 2004. 

[49] K. Knudsen and A. Tamasan. Reconstruction of less regular conductivities in the plane. Communieations 
in Partial Differential Equations, 29:361-381, 2004. 

[50] A. Lechleiter and A. Rieder. Newton regularizations for impedance tomography: convergence by local 
injectivity. Inverse Problems, 24:065009, 2008. 

[51] L. Liu. Stability Estimates for the Two-Dimensional Inverse Conduetivity Problem. PhD thesis. University 
of Rochester, 1997. 

[52] J.L. Mueller and S. Siltanen. Direct reconstructions of conductivities from boundary measurements. SIAM 
Journal on Seientifie Computing, 24(4): 1232-1266, 2003. 

[53] J.L. Mueller and S. Siltanen. Linear and Nonlinear Inverse Problems with Praetieal Applieations. SIAM, 

2012. 

[54] D. Mumford and J. Shah. Optimal approximations by piecewise smooth functions and associated varia¬ 
tional problems. Communieations on pure and applied mathematies, 42(5):577-685, 1989. 

[55] E. K. Murphy, J. L. Mueller, and J. G. Newell. Reconstructions of conductive and insulating targets using 
the D-bar method on an elliptical domain. Physiologieal Measurement, 28(7):S101-S144, 2007. 

[56] A. Nachman, J. Sylvester, and G. Uhlmann. An n-dimensional Borg-Levinson theorem. Communieations 
in Mathematieal Physies, 115:595-605, 1988. 

[57] A. 1. Nachman. Global uniqueness for a two-dimensional inverse boundary value problem. Annals of 
Mathematies, 143:71-96, 1996. 

[58] R.G. Novikov. A multidimensional inverse spectral problem for the equation + (t'(x) — eu(x))f> = 0. 
Eunetional Analysis and Its Applieations, 22(4):263-272, 1988. 

[59] G. D. Perttunen, D. R. Jones, and B. E. Stuckman. Lipschitzian optimization without the lipschitz con¬ 
stant. Journal of Optimization Theory and Applieation, 79(1): 157-181, October 1993. 

[60] Thomas Pock, Antonin Ghambolle, Daniel Gremers, and Horst Bischof. A convex relaxation approach for 
computing minimal partitions. In Computer Vision and Pattern Reeognition, pages 810-817, 2009. 

[61] Luca Rondi and Eadil Santosa. Enhanced electrical impedance tomography via the Mumford-Shah func¬ 
tional. ESAIM Control Optim. Cale. Var., 6:517-538, 2001. 

[62] B. Sandberg and T.E. Ghan. A logic framework for active contours on multi-channel images. Journal of 
Visual Communieation and Image Representation, 16(3):333-358, 2005. 



HYBRID SEGMENTATION AND D-BAR METHOD FOR EIT 


23 


[63] S. Siltanen, J. Mueller, and D. Isaacson. An implementation of the reconstruction algorithm of A. Nachman 
for the 2-D inverse conductivity problem. Inverse Problems^ 16:681-699, 2000. 

[64] J. Sylvester and G. Uhlmann. A global uniqueness theorem for an inverse boundary value problem. Annals 
of Mathematies, 125:153-169, 1987. 

[65] Cheng Tai, Xiaoqun Zhang, and Zuowei Shen. Wavelet frame based multiphase image segmentation. SIAM 
J. Imaging Sei., 6(4):25212546, August 2013. 

[66] K. van den Doel and U.M. Ascher. On level set regularization for highly ill-posed distributed parameter 
estimation problems. Journal of Computational Physies, 216(2):707 - 723, 2006. 


E-mail address: sarah.hamilton@marquette.edu 
E-mail address: reyes.juanmanuel@gmail.com 
E-mail address: samuli.siltanen@iki.fi 
E-mail address: xqzhang@sjtu.edu.cn 



